Fast computation of magnetostatic fields by Non-uniform Fast Fourier Transforms 
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The bottleneck of micromagnetic simulations is the computation of the long-ranged magnetostatic 
fields. This can be tackled on regular A^-node grids with Fast Fourier Transforms in time TV log A'', 
whereas the geometrically more versatile finite element methods (FEM) are bounded to A'''*'''^ in the 
best case. We report the implementation of a Non-uniform Fast Fourier Transform algorithm which 
brings a NlogN convergence to FEM, with no loss of accuracy in the results. 



The power of computers steadily increases over the 
years while the size of devices used in fundamental 
science or technology is shrinking. Today we have 
reached a cross-over where numerical simulations are 
capable of describing in detail the physics of nanodevices, 
for which they thus play a leading role in their unders- 
tanding and designing. In numerical micromagnetics for 
spin electronics, a bottleneck is the computation of the 
magnetostatic interactions which by nature are long- 
ranged. These interactions can be expressed in terms of 
either magnetostatic field H or scalar pseudo-potential 
<j> such that H = —Vcj). The latter is convenient since 
it boils the problem down to a single scalar unknown. 
To deal with magnetostatic interactions, essentially two 
distinct approaches have been implemented, depending 
on the type of mesh used : 



for Finite Difference (ie. translation invariant) 
meshes, a Green approach with Fast Fourier 
Transforms, called FD-FFT. Computation time is 
moderate (A^logiV with N the number of nodes), 
but the curved boundaries that occur often in 
experimental devices are not ideally described ; 

for Finite Element (much more general) meshes, 
a Finite Element Method coupled to a Boundary 
Element Method, called FEM-BEM. This can 
faithfully describe curved boundaries but compu- 
tation time is higher, at least N^/^ in 2D (rcsp. 
TV^/^ in 3D). 



In this Letter we report the implementation of a new 
magnetostatic code which combines the advantages of 
both cited approaches : it uses a FEM mesh thus des- 
cribing curved boundaries as well as FEM-BEM, albeit 
with computation time A^logiV. It is based on a algo- 
rithm reported recently for computing non-periodic Fast- 
Fourier Transforms (NFFT)i. Our code, called FEM- 
NFFT, proves to be significantly faster than FEM-BEM 



with no loss of precision. As a first step the implementa- 
tion was done for a 2D geometry (which pertains to 3D 
systems with one direction of translational invariance, i.e. 
cylinder- like) for the proof of concept. The gain is expec- 
ted to be even greater in more realistic 3D calculations. 
This demonstrates the potential of FEM-NFFT for mi- 
cromagnetism and thus spin-electronic devices. 

Let us recall the principle and features of FD-FFT and 
FEM-BEM before presenting our approach and results. 
We consider a system J7 with boundary dQ, displaying a 
known magnetization distribution M(r). 

Finite Difference micromagnetic codes use a transla- 
tion invariant grid. On such a grid, FFTs can be used to 
compute convolutions in time A^logA^. This motivates a 
Green approach for magnetostatics : is calculated as a 
convolution of the Green function G = — (l/27r)logr in 
2D [rcsp. G = l/(47rr) in 3D] with the magnetic charges, 
volumic p = — V ■ M and surfacic cr = M • n : 



(r)= /p(r')-G(r-r')dr'+ / a{r')-G{r-v')dr' (1 
Jn Jdn 
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The iVlogA^ speed explains the wide and lasting use 
of FD in micromagnetic simulations^. However, most 
devices have curved boundaries, either by design or as 
a result of experimental imperfections. Simulating these 
cases with FD requires the use of saw-tooth boundaries 
to describe the magnetic material. This geometrical 
approximation may induce inadequate descriptions'^. 

Finite Element micromagnetic codes use, on the 
contrary, complex-shaped meshes with triangles (resp. 
tetrahedrons) as 2D (resp. 3D) unit cells. They conse- 
quently suffer much less from the above-mentioned limi- 
tations. However, without translational invariance of the 
mesh, FFTs are not available. Bearing in mind that direct 
summation of Eq.(l) on a FEM mesh, called FEM-direct, 
would cost iV^ time, one understands why the Green ap- 
proach is thought incompatible with FEM codes. 

To deal with magnetostatics, FEM codes thus go back 
to the Poisson equation — A0 = p, with the usual re- 
gularity condition that the field should decay at infinity. 
Applying FEM, the equations are translated into a linear 
system, which is solved by standard iterative methods'^ 



in time N^''^ in 2D (resp. N'^/^ in 3D) Not only is this 
asymptotically slower than the TV log N time required for 
FD-FFT, but N here takes higher values, because the 
mesh must extend well beyond fi in order to tackle the 
regularity condition at infinity. This induces an additio- 
nal slowdown, and also creates finite-size artifacts. 

To avoid meshing outside fi, the main approach 
consists in coupling FEM with a Boundary Element 
Method, resuhing in the so-called FEM-BEM^. The 
asymptotic complexity of the Poisson solver is unchanged 
but the BEM step introduces another time limitation 
Nq where Nq is the number of boundary nodes. In the 
most favorable case consisting of compact systems Nq « 
7V1/2 in 2D (resp. Nq ~ iV^/^ in 3D). However for fiat 
geometries, of particular relevance to applications, Nq k, 
N, in which case the time limitation for FEM-BEM may 
be pretty severe. 

Our innovation is to revert, within the FEM frame- 
work, to a Green approach. The convolution (1) is dis- 
cretized in a way typical for FEM, and computed using 
a fairly recent mathematical method called NFFT (Non- 
uniform Fast Fourier Transform). NFFTs allows one to 
compute discrete convolutions in time N log N without 
the equispaced data requirement of FFTs. 

More in details, we seek (p at the nodes (r.^), i=i...N of 
the mesh. A linear interpolation inside each element, of 
known magnetization values M(ri), is used to evaluate 
charges p and a at points r-,-, j=i...M defined as the qua- 
drature points for the integrals in (1)^. 

Consequently, (1) is rewritten the following way : 



M 



{ri) = Y^ p,j Q{vi - Yj) uj det J{rj) 



j=i 



M 



^ cTj G (r; - Tj ) u}j det J{rj ) (2) 



where 



- pj = p{rj) if rj is in the interior of fl, otherwise, 

- aj — crirj) if rj is on dfl, otherwise, 

- tOj is the weight of tj in the quadrature scheme, 

- J{rj) is the Jacobian of the affine transformation 
mapping the unit element on that containing rj . 

We then use NFFTs to compute (2). Although the se- 
minal papcr^ dates back to 1993, the NFFT method re- 
mains little-known even in the mathematical community. 
A presentation of the method can be found Rcf.*. Here 
we sketch the basic strategy and show that computation 
time does not exceed iVlogA^, ie. that of the classical 
FFT. 

Let us look at our non equispaced data as a sum of 
Dirac functions. The goal is to find the spectrum of this 
data function. The idea is to convey initial information 



over a regular grid, so that an FFT can be used. We 
therefore choose a regular grid X and let the data dif- 
fuse to the X nodes through convolution with a Gaussian 
function (or more generally a smooth localized function). 
The Gaussian is localized in space, so we can consider 
that each piece of data diffuses only to a fixed number of 
the nearest X nodes. Computation time of this diffusion 
step is thus proportional to N. 

The question arises how to choose the period of the X 
grid. In our ease, the answer depends on how smooth the 
Green function G is. A necessary preliminary step before 
executing the NFFT is therefore to smooth G around the 
origin ; the price to pay is an afterwards correction in the 
smoothing zone. It can be shown that a grid of size p^N 
in 2D (resp. p^N in 3D) is convenient, where p is the 
degree of smoothness chosen for G. 

We then perform, according to the initial idea, a 
FFT on the X grid, in time proportional to NlogN. 
Based on the convolution theorem, what we get is the 
Fourier coefficients of the data function multiplied by 
those of the Gaussian. Therefore, we finally divide these 
numbers by the Fourier coefficients of the Gaussian to 
get the desired spectrum. The number of divisions is 
proportional to N. As a whole, the NFFT is expected 
to behave asymptotically like NlogN, as all extra steps 
behave like N. 

We implemented FEM-NFFT to 2D test cases where 
an analytical solution 0a is available, so that errors can 
be readily estimated. Interpolation and quadrature rou- 
tines are written in C++ and the NFFT package used^ 
is in C99. For G we have chosen a smoothness degree of 
2. On each test case, we provide computation times and 
error estimates for FEM-direct, the classic FEM-BEM 
and our NFFT-based method. The computed error is 
the normalized root mean square {MsL)~^{J2i=i \4'{^i)~ 
0a(r?:)|^/TV)^^^, where Ms is the saturation magnetization 
and L the system diameter set at unity, . Computations 
are done on an Intel P4 2GHz with 1GB RAM running 
Fedora 5. 

The first test case is a disk uniformly magnetized 
along the x-axis (a cylinder in 3D space). The analy- 
tical solution is (j){x) = Mg x/2. The second case is 
the so-called magic cylinder, a circular annulus of ra- 
dii i?i , i?2 where the angle between magnetization and 
the X-axis equals twice the polar angle. The name stems 
from the uniform magnetic field thus induced in the in- 
ner region. The analytical solution inside the annulus is 
(j){r,6) = Ms r cos 9 \og{r / R2) in polar coordinates. 

Tables 1 and 2 display the numerical results for the 
two cases, respectively. It can be readily seen that FEM- 
NFFT provides results very similar to FEM-direct. The 
error induced by the NFFTs is thus negligible. Compared 
to FEM-BEM, errors are comparable for non-uniformly 
magnetized systems (see Table 2) whereas for uniform 
distributions. Green approaches, to which FEM-NFFT 
belong, are more accurate by one order of magnitude 
(sec Table 1). This is because they can treate apart volu- 



mic and surfacic charge contributions. Concerning com- 
putational time, FEM-direct is as expected the quickest, 
gaining a factor around 5 over FEM-BEM for the finer 
meshes. Nodes required in 3D cases of interest commonly 
count up to 10^, around which number the time advan- 
tage of FEM-NFFT over classical methods such as FEM- 
BEM is expected to reach one order of magnitude. 




Fig. 1: Test case 1 : the uniformly magnetized disk. Mesh used 
for A*' — 400 (left) and magnetization distribution (right). 



Tab. I: Test case 1 : error (in ppm) and computation time 
(in seconds) of FEM-direct, FEM-BEM and FEM-NFFT for 
different mesh sizes. 



iV 


error 


time 1 


FEM- 


FEM- 


FEM- 


FEM- 


FEM- 


FEM 




dircct 


BEM 


NFFT 


direct 


BEM 


NFFT 


400 


93.6 


162 


90.3 


0.401 


0.070 


0.089 


1572 


17.3 


52.8 


16.8 


7.71 


0.390 


0.385 


9489 


2.16 


25.8 


2.09 


233 


4.75 


2.18 


37938 


— 


23.2 


0.650 


3720* 


40.1 


9.81 




To conclude, we have successfully implemented a Non- 
uniform Fast Fourier Transform (NFFT) algorithm to 
compute magnetostatic fields for micromagnctic simula- 
tions based on Finite Element methods (FEM) . The new 
approach, called FEM-NFFT, combines the advantages 
previously found separately in Finite Difference methods 
(computation time scaling like A^log A^) and FEM (faith- 
ful description of curved boundaries). Thus FEM-NFFT 

Tab. II: Test case 2 : error (in ppm) and computation time 
(in seconds) of FEM-direct, FEM-BEM and FEM-NFFT for 
different mesh sizes. 



iV 


error 


time 1 


FEM- 


FEM- 


FEM- 


FEM- 


FEM- 


FEM 




direct 


BEM 


NFFT 


direct 


BEM 


NFFT 


1192 


174 


172 


174 


3.58 


0.290 


0.331 


2091 


95.8 


94.2 


95.8 


14.0 


0.640 


0.498 


8214 


24.7 


32.1 


24.8 


211 


4.59 


2.02 


22683 


— 


21.5 


9.05 


1600* 


34.0 


6.62 



estimated. 



Fig. 2: Test case 2 : the magic cylinder. Mesh used for TV 
1192 (left) and magnetization distribution (right). 



promises a leap in the attractiveness of micromagnetic 
simulations of spin electronic devices. 
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